#include <iostream>
#include <fstream>
#include "MultigridSolver.h"
#include <cmath>
#include <algorithm>
#include <vector>
#include <jsoncpp/json/json.h>

#define func F1

template <int dim>
class F1 : public Function<dim>
{
public:
    double operator()(double x) override
    {
        return exp(sin(x));
    }
    double diff(double x) override
    {
        return exp(sin(x))*cos(x);
    }
    double Laplace(double x) override
    {
        return exp(sin(x))*sin(x)-exp(sin(x))*cos(x)*cos(x);
    }
};

template <int dim>
class F2 : public Function<dim>
{
public:
    double operator()(double x) override
    {
        return -exp(-x);
    }
    double diff(double x) override
    {
        return exp(-x);
    }
    double Laplace(double x) override
    {
        return exp(-x);
    }
};

template <int dim>
class F3 : public Function<dim>
{
public:
    double operator()(double x) override
    {
        return cos(x)+1;
    }
    double diff(double x) override
    {
        return -sin(x);
    }
    double Laplace(double x) override
    {
        return cos(x);
    }
};

int main()
{
    Json::Reader reader;
    Json::Value root;
    std::ifstream in("hw04.json", std::ios::binary);
    if (!in.is_open())
      {
		std::cout << "Error: cannot open the file!" << std::endl;
		return -1;
      }

    int maximum_iteration_number;
    double relative_accuracy;
    std::vector<int> n;
    std::string boundary_condition, restriction_operator, interpolation_operator, grid_type, initial_guess;

    if (reader.parse(in, root))
      {
        for (int i = 0; i < root["n"].size(); i++)
	    {
	        int temp_n = root["n"][i].asInt();
	        n.push_back(temp_n);
	    }
        boundary_condition = root["Boundary_condition"].asString();
        restriction_operator = root["Restriction_operator"].asString();
        interpolation_operator = root["Interpolation_operator"].asString();
        grid_type = root["Grid_type"].asString();
        initial_guess = root["Initial_guess"].asString();
        maximum_iteration_number = root["Maximum_iteration_number"].asInt();
        relative_accuracy = root["Relative_accuracy"].asDouble();
      }

    const int dim = 1;
    func<dim> f;
    std::vector<double> Error_norm;
    for (int i = 0; i < n.size(); i++)
      {
        std::vector<double> v;
        if (dim == 1 && initial_guess == "Zero") 
        {
            v = std::vector<double>(n[i]+1, 0.0);
        } 
        Solver<dim> V(f,boundary_condition, restriction_operator, interpolation_operator, maximum_iteration_number, relative_accuracy, 1.0/n[i], v);
        if (grid_type == "VC")
        {
            V.Vcycle();
        }
        if (grid_type == "FMG")
        {
            V.FMG();
        }

        std::vector<double> R = V.get_residual();
        double error_norm = V.get_errornorm();
        Error_norm.push_back(error_norm);
        R.erase(R.begin());

        std::cout << "n = " << n[i] << std::endl;
        std::cout << "Iteration number: " << R.size() << std::endl;
        std::cout << "Residual's norm & Relative reduction rate: " << std::endl;
        std::cout << R[0] << std::endl;
        for (int j = 1; j < R.size(); ++j)
        {
            std::cout << R[j] << " & " << round(R[j]/R[j-1]*1000)/1000 << std::endl;
        }
      }

    std::cout << "Errornorm of each n: " << std::endl;
    for (int i = 0; i < Error_norm.size(); ++i)
    {
        std::cout << Error_norm[i] << ", ";
    }
    std::cout << std::endl;
    std::cout << "Relative convergence rate(1og 2): " << std::endl;
    for (int i = 1; i < Error_norm.size(); ++i)
    {
        std::cout << log2(Error_norm[i-1]/Error_norm[i]) << ", ";
    }
    std::cout << std::endl;

    return 0;
}
